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CN ' Abstract. 

HST observations reveal that young massive star clusters form in gas-rich environments like the 
Antennae galaxy which will merge in collisional processes to form larger structures. These clusters 
amalgamate and if some of these clusters harbour a massive black hole in their centres, they can 
become a strong source of gravitational waves when they coalesce. In order to understand the 
dynamical processes that are into play in such a scenario, one has to carefully study the evolution of 
the merger of two of such young massive star clusters and more specifically their respective massive 
black holes. This will be a promising source of gravitational waves for both, LISA and the proposed 
Big Bang Observer (BBO), whose first purpose is to search for an inflation-generated gravitational 
Q\ \ waves background in the frequency range of 10~ 1 — 1 Hz. We used high-resolution direct summation 

t • . N— body simulations to study the orbital evolution of two colliding globular clusters with different 

initial conditions. Even if the final eccentricity is almost negligible when entering the bandwidth, it 
will suffice to provide us with detailed information about these astrophysical events. 
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Nowadays it is well established that massive stellar clusters form in interacting galax- 
ies. High-resolution Hubble Space Telescope Observations of the Antennae (J, 0] or 
Stephan's Quintet [3] show that hundreds of young massive star clusters are lurking in 
the star forming regions and that they are clustered into larger clusters of a ~ few 100 pc. 
The images also reveal that only a few of these clusters are reddened, thus suggesting that 
the gas has been already removed in most of them. Since they harbour ~ 10 5 stars within 
~ 1 — few parsecs and are older than 5 Myr, they are most likely to be bound clusters. Su- 
pernova explosions originating from stars with a mass larger than 8M only contribute 
to the total cluster mass -which for ~ 10 7 — 10 5 stars is about ^# c i ~ 10 7 — 10 5 M Q - in 
about ~ 10%. Stars more massive than 2M & contribute in only ~ 25% to the total, so 
that the clusters are very likely to be bound. 

The so-called clusters complexes have been employed in the literature as a possibility 
to build compact dwarf galaxies as a result of the amalgamation of their smaller clusters 
constituents fl, @] in collisional processes. On the other hand, it has been studied by 
different authors how in such a young dense cluster collisional processes among the 
heaviest stars that segregate to the centre because of dynamical friction might lead 
to the formation of a very massive star [@, 0, HI @]. Such star will become unstable 
and thus possibly create a massive black hole in their centre with a mass ranging 
between 10 2 — 1O 4 M -which therefore receives the surname 'intermediate-mass' black 
hole (IMBH)-. This means that the young dense clusters in the cluster complexes are 



FIGURE 1. Initial set up for a parabolic collision in the centre-of-mass of both clusters 



possibly guarding an IMBH in their centres; since the clusters merge with each other, 
the respective IMBHs will have, at least, the possibility of forming a bound system. 
Whether or not they will merge within a Hubble time and what the implications for 
LISA and the BBO would be, is something to be analysed numerically. 



HOW TO MERGE TWO YOUNG DENSE CLUSTERS 

The clusters were set on to a parabolic orbit so that the minimum distance at which they 
pass by is of Fig. © if they are considered to be a point particle at that moment. In 
the centre of mass reference frame of both clusters, as shown in the figure, we have that 
xi = A^d, X2 = —X\ d, vi = A2 v re i and V2 = — Ai v re i; where v re i is the relative velocity 
of the clusters, Xi 2 their positions (if we regard them to be a point mass, or to their 
centres) and Ai 5 2 = (mi +^2) 

From the reduced particle standpoint, we have to determine d and v re i when the 
separation is d (which is given as an initial condition). Since the reduced mass is 
fl = m\ m i/ {mi + mj), for a parabolic orbit we have that the energy at the pericentre 

-Gm x m 2 1 2 ~Gm\ m 2 1 , ... 
E = —j— + ^ = — ^ + -M vLx = 0. (1) 

Since the specific angular momentum per unit /1 is I = l z = — v max ■ d m i n = |d A v re i| z = 
xv y — yv x , we obtain v max = 

\/2G{m\ +m 2 )/d min , for a given specific angular momen- 
tum I = ^JlGd m i n {m\ +ni2). 

As for the components of the velocities, for a given relative velocity of v re i = 
A/2GJ m i n (mi +m.2)/d and with the help of the relations — / = x ■ v y — y ■ v x , v^ el = v 2 x + Vy 
and d 2 = x 2 +y 2 , we can infer that the velocity components are 



We resort finally to the definition of parabola to obtain the required expressions for x 
and y, 2 d m \ n — x = d,so that x = d — 2d m i n and y = V d 2 —x 2 . 

For the model presented in this work, we chose a system with d m [ n = 2 pc, correspond- 
ing to a relative velocity at pericentre of 23.3 km/s. The clusters will always merge and 
form a larger cluster, because they are initially set in a parabolic orbit. The number of 
stars used for each cluster is Jf* = 6.3 x 10 4 , the masses of they set to 6.3 x 1O 4 M and 
we used for the initial distribution a King model of concentration Wq = 1 (cluster 1) and 
Wq = 6 (cluster 2). The central velocity dispersions were set to <T CO rel =8.41 km/s and 
0core2 = 8.29 km/s (hereafter the subscripts 1 and 2 stand for the cluster 1 and 2), with 
core radii of R COK \ = 0.203 pc and R CO k2 = 0.293 pc. Both clusters host additionally an 
EVIBH of mass 300M Q in their centres. 

The simulation was performed on a special-purpose hardware GRAPE-6A single PCI 
card with a peak performance of 130 Gflops 1 1011 - roughly equivalent to 100 single 
PCs, with the direct-summation NBODY4 code of Aarseth 111 ill . This choice was made 
for the sake of the accuracy of the study of the orbital parameters evolution of the 
binary of IMBHs; for this numerical tool includes both the KS regularisation and chain 
regularisation, which means that when two or more particles are tightly bound to each 
other or the separation among them becomes very small during a hyperbolic encounter, 
the system becomes a candidate to be regularised in order to avoid problematical small 
individual time steps. The basis of direct NBODY4 codes relies on an improved Hermit 
integrator scheme lfl2ll for which we need not only the accelerations but also their 
time derivative. The computational effort translates into accuracy and this way we can 
reliably follow track of the orbital evolution of every single particle in our system. Other 
alternative codes that add a softening to the gravitational forces (i.e. substituting the l/r 2 
factor with l/(r 2 + e 2 ), where r is the radius and £ the softening parameter) in order to 
avoid them to become too large are to be discarded if we want to befittingly make a 
highly accurate estimate of the orbital evolution of the IMBHs system (for at a certain 
point in the evolution of the binary the separation ~ £) which is the final purpose of our 
numerical study. 



DYNAMICS AND GEOMETRY OF THE SYSTEM 

In Fig. © we show nine snapshots of the merger of the two clusters and the position 
of the IMBHs marked with a black dot. The nine snapshots show the evolution of the 
clusters merger for T = 0, 54, 66, 76, 110, 128, 146, 152 and 192 ^tnb N-body time 
units, which correspondingly are T = 0,2.17, 3.19, 4.62, 5.37, 6.13, 6.38 and 8.06 
Myrs. We can clearly observe an exchange of stars between the clusters already after 
the first cluster interaction (the corresponding "red" stars do not show up in the "blue" 
cluster because the blue colour overwrites the red). As shown in Fig. ©, the binary 
of intermediate-mass black holes hardens steadily. We obtain the classical value of the 
hardening constant H ~ 16 IU3LI14I1 . 

In Fig. © we show the evolution of the triaxiality of the cluster formed as a result 
of the merger of the two cluster for our fiducial model. In the figure, a, b and c -with 
a > b > c ab definitio- are the semi-major axes of the ellipsoid of inertia, determined by 
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FIGURE 2. Projection in the X-Y plan of the trajectory of every single particle in the process of merger 
of the two young dense clusters. See text for further details 
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FIGURE 3. Evolution of the semi-major axis of the massive black hole binary 



four different mass fractions of the stars which were distributed according to the amount 
of gravitational energy. This means that the lower the mass fraction is, the closer we are 
to the centre of the resulted merged system. We can clearly see that the system resulting 
of a realistic parabolic cluster merger has a shorter axis and two approximately equal 
longer axes, but larger than the third one, a ~ b > c, which means that it is manifestly 
oblate and not prolate. We can also observe that the binary of IMBHs makes the system 
more spherical in the centre than at larger distances, where it is flatter. This is so because 
the flattening depends on the gravitational and centrifugal forces of the cluster, its size, 




FIGURE 4. Triaxiality of the resulting merged cluster for different mass fractions (left panel) and for 
the mass fraction 0.5 (right panel) 



rotation and density of stars. For the case of two galactic nuclei harbouring massive black 
holes, this would translate into a clear danger of hang-up of the binary or stalling; i.e. 
there would not be enough stars to strongly interact with the binary of massive black 
holes because the loss-cone -a phase space region where stars can interact with the 
central object/s in one crossing time, II15L 11711 - would remain depleted because of 
the lack of centrophilic orbits in an axially symmetric potential. In the case studied here, 
though, this is not the case, because the ratio of relaxational to dynamical timescales is 
small due to the lower number of stars in the system. 



IMPLICATIONS FOR LISA AND THE BBO: GENERAL 
DISCUSSION ON NBODY4 PROSPECTS FOR 
LOWER-FREQUENCY GRAVITATIONAL-WAVES 
ASTROPHYSICS 

After a strong dynamical evolution with the surrounding stars of the merged cluster 
in which the IMBHs are embedded, the evolution is dominated by gravitational wave 
emission during the last 10 8 yrs. The evolution from the moment the binary is bound 
takes ~ 160 Myrs and enters the LISA bandwidth (i.e. has an orbital period of less than 
10 4 sec) on an almost circular orbit (for a more detailed description of the final stage of 
the evolution see [18]). For the work presented here we performed also a parallel test 
simulation including the possibility that stars are tidally disrupted by the IMBHs binary 
if they enter the tidal disruption radius of one of the IMBHs but these effects are almost 
insignificant in terms of the description of the parameter space of the IMBH binary and 
thus negligible for the main purpose of this study. In any case, the "final" eccentricity 
of the binary (the eccentricity it achieves by entering the LISA band) is almost zero; 
i.e. the binary has almost totally circularised. We can expect an event rate of between 



four and five of these events per year for LISA lfl8ll . As for the detectability of the 
gravitational waves emitted by the IMBH binary as it inspirals and merges, in H 1 811 we 
show the amplitude ratio assuming a coherent integration (matched filtering) over the 
observation time. For both observatories, LISA and the BBO, only the n = 2 harmonic 
of the quadrupolar emission is clearly detectable during the last few years of inspiral at 
1 Gpc, since it is dominant for low eccentricity binaries; some higher (odd) harmonics 
depend on the mass ratio and are zero for equal mass (for zero eccentricity orbits). The 
eccentricity is so low (10~ 3 ) that the source should be closer than ~ 50— 100 Mpc to be 
detectable. 

For the data analysis of these sources, it is very important to know the phase of the 
gravitational wave; a missmatch even in one single cycle can reduce the signal to noise 
ratio (SNR). In practice, since we do not know the parameters of the gravitational wave 
signal, we will search over a large parameter space (including the masses of the objects, 
the eccentricity etc) and the template with largest SNR will be chosen. The highest SNR 
will be achieved by resorting to templates with parameters close to the actual ones. 
However, there are some correlations between the parameters, so that one can compute 
the level of correlation and the error bar on the parameter determination. For strong 
signals with long duration within the detector band, we need to know the waveform 
(more specifically the phase) very accurately. The eccentricity affects the phase as well 
as the amplitude. Since the analysis is very sensitive to even a tiny change in phase, we 
can expect that one will be able to subtract the necessary information of the waveforms 
of the event to identify it, due to the induced detectable phase shift of the residual 
eccentricity (Bernard Schutz personal communication). 

The space-borne LISA mission will fly in ~ 10 yrs and critical design choices affect- 
ing the ability to detect this kind of events will have to be made soon. It is of big im- 
portance to produce robust estimates for the rates and typical orbital parameters of these 
and other events interesting for the detection of gravitational waves (extreme mass-ratio 
inspirals, for instance) in order to develop a detection template family. This detection 
template family consists of a bank of waveforms. It is relatively fast to generate wave- 
forms depending on the initial conditions and mass ratio, but not fast enough if one 
wants to have a bank of them. If we want to rigorously explore the parameter space of 
these events, we need realistic Astrophysical estimates of the eccentricity, mass ratio etc 
at the beginning of the final merger, "one step before" the objects enter the LISA band. 
An assumption for the initial parameter space is necessary in order to develop wave- 
form banks. For instance, in the case of an inspiral search, we know that majority -if 
not all- of binaries entering the bandwidth of ground-based detectors like LIGO will be 
circular. This allows us to significantly simplify the study of the parameter space. Sim- 
ilarly, in the case of space-born detectors like LISA, we need to better understand what 
we will observe and so focus our efforts on the part of the parameter space that from 
which we expect the largest contribution, even if our final aim is to cover the whole pa- 
rameter space. For this, an Astrophysical understanding of the scenario is of paramount 
importance. The high-accuracy, Astrophysical numerical simulations performed with a 
direct- summation NBODY4 will shed light on many of these aspects by studying realis- 
tic astrophysical scenarios with relativistic corrections, which have been added recently, 
and arbitrary geometries 1 1^, 11]. The biggest impact on data analysis is a realistic esti- 
mation of the event rate above a fixed SNR and possibly also as a function of the SNR 



itself. The models one can develop with these numerical tools have no precedents be- 
cause of the inclusion of non-symmetry, rotation and relativistic corrections. A realistic 
estimation of the parameter space is crucial for the data analysis of LISA gravitational 
waves. A too small event rate and signals with smaller SNR will require a method sim- 
ilar to ground based gravitational waves analysis. This will imply the need for digging 
out of signals from the noise with only a few overlapping signals in the data stream. On 
the other hand, in the case of having a large even rate, many overlapped signals would 
probably require different data analysis algorithms to make a parameter estimation. 



We acknowledge the Astronomisches Rechen-Institut for the computing resources on 
the GRACE cluster of the Volkswagen Foundation, SFB439. The work of PAS has been 
supported in the framework of the Third Level Agreement between the DFG and the IAC 
(Institute de Astroffsica de Canarias). The author thank the Kavli Institute of Theoretical 
Physics for inviting him to attend the Physics of Galactic Nuclei program, where he 
finished the manuscript. The work was thus partially supported by the NCF under Grant 
PHY99-07949. He also shows gratitude to Stas Babak, Bernard Schutz and Curt Cutler 
for enlightening conversations. The work reported here has been done in collaboration 
with Marc Freitag (IoA, Cambridge) 



1 . B.C. Whitmore, Q. Zhang, C. Leitherer, S. M. Fall, F. Schweizer, and B. W. Miller, ApJ 118, 155 1- 
1576 (1999), |asFro-ph/9907430l 

2. Q. Zhang, and S. M. Fall, ApJL 527, L81-L84 (1999), astro-ph/9911229| 

3. S. C. Gallagher, J. C. Charlton, S. D. Hunsberger, D. Zaritsky, and B. C. Whitmore, ApJ 122, 163- 
181 (2001), |asFro-ph/ 01 040051 

4. M. Fellhauer, and P. Kroupa, MNRAS 330, 642-650 (2002), [ast ro-ph/ 0110 621 1 

5. M. Fellhauer, and P. Kroupa, MNRAS 359, 223-227 (2005), |astro-ph/0502 08 6| 

6. S. F. Portegies Zwart, and S. L. W. McMillan, ApJL 528, L17-L20 (2000), |astro-ph/991006Tl 

7. M. A. Giirkan, M. Freitag, and F. A. Rasio, ApJ 604, 632-652 (2004). 

8. S. F. Portegies Zwart, H. Baumgardt, P. Hut, J. Makino, and S. L. W. McMillan, Nature 428, 724-726 
(2004), [ast ro-ph/ 4 2 6 2 2| 

9. M. Freitag, M. A. Giirkan, and F. A. Rasio, MNRAS 368, 141-161 (2006), |ast r o-ph / 5 3 1 3 0| 

10. T. Fukushige, J. Makino, and A. Kawai, ASP 57, 1009-1021 (2005). 

11. S. J. Aarseth, Gravitational N-Body Simulations, ISBN 0521432723. Cambridge, UK: Cambridge 
University Press, November 2003., 2003. 

12. S. J. Aarseth, ASP 111, 1333-1346 (1999). 

13. G. D. Quinlan, NewA 1, 255-270 (1996). 

14. A. Sesana, F. Haardt, P. Madau, and M. Volonteri, ApJ 611, 623-632 (2004), 
|astro-ph/ 04 015431 

15. J. Frank, and M. J. Rees, MNRAS 176, 633-647 (1976). 

16. P. Amaro-Seoane, and R. Spurzem,A«V/M5 327, 995-1003 (2001). 

17. P. Amaro-Seoane, M. Freitag, and R. Spurzem, MNRAS 352, 655-672 (2004), 

|astro-ph/0401163| " 

18. P. Amaro-Seoane, and M. Freitag, Accepted for publication ApJL (2006), ast ro-ph /0610478| 

19. G. Kupi, P. Amaro-Seoane, and R. Spmzem,MNRAS pp. L77+ (2006), [ast ro-ph / 6 2 1 2 5| 



ACKNOWLEDGMENTS 



REFERENCES 



